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Abstract 

We investigate superradiance of an ensemble of atoms coupled to an integrated superconduct- 
ing LC-circuit. Particular attention is paid to the effect of inhomogeneous coupling constants. 
Combining perturbation theory in the inhomogeneity and numerical simulations we show that 
inhomogeneous coupling constants can significantly affect the superradiant relaxation process. In- 
complete relaxation terminating in "dark states" can occur, from which the only escape is through 
individual spontaneous emission on a much longer time scale. The relaxation dynamics can be sig- 
nificantly accelerated or retarded, depending on the distribution of the coupling constants. On the 
technical side, we also generalize the previously known propagator of superradiance for identical 
couplings in the completely symmetric sector to the full exponentially large Hilbert space. 



1 



I. INTRODUCTION 



With the advance of experimental quantum information processing, the need to study 
"hybrid quantum processors" has arisen l|-|8|]. In such a processor different physical systems 
are used in order to exploit their respective advantages, such as long coherence times versus 
fast processing times, or fast propagation in the case of quantum communication. A natural 
candidate for a hybrid quantum processor is a system of atoms coupled to superconducting 
circuits [3]. Circuit-QED schemes are rapidly emerging as a promising new avenue towards 
scalable quantum computation. These schemes combine the strong coupling and precise 
control of cavity-QED systems with the scalability of integrated solid state circuits. Nev- 
ertheless, the coherence times achieved so far are of the order of a /is for transmon qubits 
9] . Much longer coherence times are achievable for superpositions of the hyperfme levels of 
alkali-metal atoms, a fact well known from the design of atomic clocks. Therefore, atoms are 
predestined as a building block for a quantum memory. The coupling of the hyperfme levels 
to their environment is through a magnetic dipole transition and thus much weaker than the 
coupling through an electric dipole. This limits the bandwidth with which information can 
be transferred to and from such a quantum memory. It is beneficial to use a large number of 
atoms to store a single excitation as the single photon Rabi frequency that determines 
the rate of energy exchange scales as \/N with the number of atoms N. 

Beyond their technical relevance, hybrid quantum processors are also interesting for 
studying physical effects on a more fundamental level. They allow for easily modifiable 
parameters, and to reach new parameter regimes, so far inaccessible in traditional quan- 
tum optical systems. For example, it has been proposed that the ultra-strong couplings 
of circuit-QED schemes might allow one to observe the phase transition in superradiance 
predicted last century by Mallory, and Hepp and Lieb 1CH12|. The existence of such a 



transition was subject of considerable theoretical debate (see e.g. Il3|), culminating in the 
recent claim that the transition can in principle not be observed in cavity-QED, but should 
be observable in circuit-QED [jjl . For atoms coupled to a superconducting circuit one might 
envisage to perform precisely controlled experiments on superradiance. Superradiance is a 
rather complex effect that can include things like mode competitions or beating 15|. It 
would therefore be desirable to go beyond previous experiments 
cavity QED in terms of control of different parameters. 



15Ml9| using "traditional" 
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FIG. 1: (Color online) N two level atoms in a dipole trap with a transition in resonance with a 
superconducting LC-circuit. 

In this paper we study superradiance of an ensemble of atoms coupled to an integrated 
resonant LC circuit. We show that there is indeed a regime in which superradiant behavior is 
expected. At the same time, additional complications arise due to variations in the coupling 
constants, which are fundamentally related to the small size of the integrated LC circuit. 
Most of this paper is therefore dedicated to the study of the effects of inhomogeneous coupling 
constants on superradiance. So far the study of superradiance has been almost entirely 
limited to homogeneous coupling constants, and/or initial states that are fully symmetric 
under all permutations of atoms. We develop a general theoretical framework that allows 
one to deal with such inhomogeneities. Based on perturbation theory in the inhomogeneity, 
this framework starts with the construction of a semiclassical solution of the superradiant 
master equation in the entire exponentially large Hilbert space for the homogeneous case. 



II. MODEL 



A. Physical system 

We consider N atoms that are, on the time scale of the experiment, held at fixed positions 
close to an LC-circuit (see Fig. [TJ). We assume that a single atomic transition is in resonance 
with a single mode of the cavity. To be specific, we are particularly interested in an ensemble 
of 87 Rb atoms trapped and cooled in a dipole trap in close vicinity (a few /im) of the 
surface of the LC-circuit. Neutral 87 Rb has a hyperfine split ground state with total angular 



momentum F = 1 and F = 2. The hyperfine sp_ 



1 and F = 2 is u ~ 2vr x 6.834 GHz 



20 



itting between the states \F,m F ) with 
2l\ . Without an external magnetic field, 



the Zeeman sublevels with ^-component of the total angular momentum mp = —F, . . . , F 
are degenerate. We will focus on the situation where an additional small static B-field 
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in the z— direction is applied that splits the degenerate hyperfine states, and defines the 
quantization axis for the atoms. Neglecting the magnetic dipole moment of the nucleus, and 
for small magnetic field strength B <^ tvjj/ (ps/xb), the magnetic moment of an atom is given 



where g s — 2 and S are the ^-factor and total electron angular momentum vector, respec- 
tively, and hb = eh/{2m e ) is the Bohr magneton (with e and m e the electron charge and 
mass) . 

We consider a simple current carrying loop with inductance L, that is part of a resonant 
LC circuit with frequency oj. The resonator can be treated as a harmonic oscillator with 
Hamiltonian Hlc = fioj(a^a + 1/2). The creation and annihilation operators and a are 
related to the current I in the circuit by / = ^/ hu / (2L)(a^ + a). That current gives rise to 
a magnetic field 



where b(x) is a dimensionless mode function with components b x , b y , b z depending on the 
geometry of the LC circuit, d a typical linear dimension of the LC circuit, and /i the 
magnetic constant. The magnetic moment of a single atom at position x couples to the 
magnetic field through the interaction Hamiltonian = — /x • B(x), which can be written 
in the |F, mp) basis as 



Starting in \F, mp) = 1 1 , 0) , the only non- vanishing resonant transition due to a B^(x) 
from the circuit oscillating at angular frequency uj is then to \F',m' F ) = |2,0) with matrix 
element (1, 0|b • S|2, 0) = hb z /2. Note that to first order in the additional static B-field, 
1 1,0) and |2,0) are not shifted by the field. In the interaction picture with respect to the 
free Hamiltonian H = H LC + fku(\2, 0) (2, 0| - |1, 0)(1, 0|)/2, the operators a and |2,0)(1,0| 
acquire time-dependent phase factors e~ lu>t and e +twt , respectively, which justifies the use 
of a rotating wave approximation. The interaction Hamiltonian for N atoms at position Xj 
then takes the familiar form 



by 



(1) 



B LC (x) = no-— b(x) , 



(2) 




N 




(4) 



i=i 
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FIG. 2: (Color online) (a) Coupling constant —g/(2Tt) as a function of the position (x,d/2,z) 
for a Rb atom close to a square loop LC resonator in the xy-plane with inductance L =lpH 
and resonance frequency uj ~ 2ir x 6.834 GHz for the transition \F,mp) = 1 1 , 0) —> |2,0). (b) 
Distribution of dimensionless coupling constants g/G with G = ^f^VME ^ or a t° ms held in 
a harmonic trap with frequencies £l x = £l y = £l z = 2ir x 1kHz centered at (d/2,d/2,d) with 
d = 10/xm, and otherwise same parameters as in (a), and temperatures T = 0.1 fiK (blue/left 
curve), T = 1/iK (black/middle curve), and T = 10 fiK (red/ right curve). 

with Qi = p(xi), 

and c_ =j |1, 0)(2, 0|j is the Pauli lowering operator for atom z. 

In Fig. [2]we show — p»/ (27r) for an LC circuit in the form of a square loop of size d = 10/mi 
in the x|/-plane, and a 87 Rb atom located at position (x, y = d/2, z) on the symmetry axis 
perpendicular to the loop, where the origin of the coordinate system is located at the lower 
left corner of the square, and the center of the square is (d/2, d/2, 0). The coupling for 
an atom at (d/2, d/2, r) with r = 5 /xm is — ~ 2tt x 388 Hz. For r = 4/xm, it increases 
to — Qi ~ 27r x 504 Hz, whereas at r = 6 /xm, we have — ^ ~ 27r x 297 Hz. We consider 
a thermal cloud of atoms confined by a 3D harmonic trap with trapping frequencies Q\ 
(A G {x,y,z}), centered at (xo,y z ,zo) with atom positions Xj taken as classical variables, 
distributed according to the probability density p(x) = p x (x — Xo)p y (y — yo)p z (z — z ). Here, 

PA(0 = ^^ex P (-e 2 AlT) (6) 
is the thermal equilibrium distribution at temperature T for a one-dimensional harmonic 
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oscillator of frequency £l\ and atomic mass M, with l\ y T = l\coth 1/2 (hQ x /(2k B T)), I 



yjh/ (MQ\), and A G {x,y,z} (see e.g. problem 2.6 in 22]). From p(x) we obtain the 
distribution of coupling constants 

P(g) = J dxdydz5(g - g>(x))p(x) , (7) 

where the integration is over the entire space R 3 . Examples of p(g) are shown in Fig. [21 
We see that even at temperatures as low as 0.1/iK, relatively large trapping frequencies of 
27rxl kHz in all directions, and a loop size of d = 10/im, the coupling constants can easily 
vary by 50% or more. For higher temperatures T = lfiK, long tails in p(g) with substantially 
larger absolute coupling constants appear due to atoms close to the loop. For T = lOyuK, the 
spread in position becomes comparable to d, and many atoms are therefore located "outside" 
the loop, leading to coupling constants with opposite sign, and a p(g) distributed around an 
average value close to zero. It is therefore clear that the consequences of the inhomogeneity 
of the coupling constants need to be investigated. 

B. Superradiance 

Superradiance is a collective emission process of iV two-level atoms resonantly coupled 
to a single mode of a resonator (under suitable conditions it can even be observed without a 



resonator, see [15j). It is observed under several conditions [15|, |23j]: 1.) The temperature of 
the electromagnetic environment must be much smaller than the level spacing of the atoms, 
such that photons only leave the resonator, whereas the entering of thermal photons can be 
neglected, and 2.) the resonator must be rather leaky. More precisely, if all N atoms couple 
identically to the resonator, one should have T <C g <C k, where k is the single photon 
escape rate from the resonator, and T denotes the single atom spontaneous emission rate. 
The latter is, in our case of hyperfme levels of 87 Rb atoms, entirely negligible. Photon escape 
from the resonator is due to the finite quality factor, Q, of the LC circuit. Q factors of up 
to 5 • 10 5 have been achieved for superconducting LC circuits [24j, leading to k = co/Q of the 
order of several kHz. For such a high-quality resonator, and the typical coupling constants 
of a few 100 Hz calculated in section III Al the number of atoms would be restricted to 
about 100. However, the quality factor of the LC circuit can be easily decreased (e.g. by 
increasing the temperature, or using a more lossy substrate), or the distance of the atoms 
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from the circuit can be increased, thus reducing the coupling constants. One can therefore 
accommodate much larger numbers of atoms. Below we also derive a lower bound on the 
number of atoms for the validity of the present analysis. 

Under the above conditions, and the assumption of weak inhomogeneity g$ = g + 5gi with 
\$9i\ ^ 9i the superradiant master equation in the rotating frame 

p = j([J-, P J+}+h.c.) (8) 

can be derived. Here, p represents the reduced density matrix describing the hyperfme 
states of the atoms, after tracing out the mode of the resonator and its electromagnetic 
environment. As mentioned above, the external state of the atoms is taken as a classical 
degree of freedom, uncorrelated from the internal states of the atoms, and with the atoms 
at fixed positions on the time scale of the experiment. The rate 7 is linked to g and k by 
7 = g 2 /n. The J± are collective ladder operators, related to the single-atom Pauli matrices 
in the Hilbert space of the two levels participating in the resonant transition by 

J± = 5><7®, (9) 

where = gj g = 1 + bg^j g denotes a dimensionless coupling strength, and g is a reference 
coupling strength (which might be chosen e.g. as the average coupling constant over all 
atoms). For later use we also define 8g~i = Sgi/g. The derivation of Eq. ([8]) follows closely 



;he calculation in 



23]. One checks that in the case of weak asymmetry the derivation in 
23] remains valid if the pseudo-angular momentum operators J± = £V a± used in [23]] are 
replaced by J± defined by Eq. fl9]). 

The homogeneous situation, with identical coupling constants, has the tremendous ad- 
vantage of the dynamics remaining restricted to a single irreducible representation (irrep) 
of SU(2) with angular momentum quantum number j. In particular, if one starts with 
all atoms excited, or any other state fully symmetric under permutation of atoms, then 
j = jmax = N/2 is the maximum possible (pseudo-) angular momentum number. Thus, 
instead of having to track 2 N states, we only have to deal with 2j max + 1 = N + 1 states. 
In the following we will consider the situation of weak inhomogeneity and assess the effect 
of the deviations 5gi to first order perturbation theory. We also introduce a rescaled dimen- 
sionless time, t = 2J , yt, J = j m ax + 1/2? which means that everything will be expressed in 



terms of the classical time scale of the system 



25| . In particular, probabilities are known to 



propagate on a time scale r ~ 1 , whereas macroscopic coherences in superradiance normally 
decay on the much shorter time scale r ~ 1/ J. 

Before studying the superradiance dynamics, let us determine the parameter regime in 
which we expect our approximation of atoms with fixed classical positions to be valid. The 
one-atom reduced density matrix of a non-interacting atom gas in thermal equilibrium in a 
ID harmonic oscillator in position representation is given by 



We see that the coherences deca y o n a length scale I c t, which is given by l Cy T 



l x tanh 1 ^ 2 (/3^r2 z /2) and /3 = Xjk^T [221 ] . Quantum effects in the external degree of freedom 
of the atoms can be neglected if Z C) t "C l x , the harmonic oscillator length scale of confinement 
of the atoms (and correspondingly for directions y,z). This implies hQ x /2 <C &bT, which 
leads to roughly 10 _2 (z/ x /kHz)/(T//iK) <C 1. Thus, for trapping frequencies of the order of 
1 kHz treating the atom positions classically is reasonable down to temperatures of order 
100 nK, and for trapping frequencies of order 100 Hz even 10 nK still means that the atom 
positions are essentially classical variables. 

A second restriction arises from the requirement of fixed atom positions on the time 
scale of the experiment. The typical thermal velocity in the x-direction of an atom follows 



from the equipartition theorem, v x = ■sJk^TjM . A typical time scale of an experiment 
is given by the decay time of the superradiant pulse, t ~ If(N'-f). During that time the 
atom should move at most a distance much smaller than d, the typical length-scale on which 
g(x) changes. This leads to a lower bound on the number of atoms, N ^> y/k^TjM / (70?). 
Together with the upper bound iV <C {^/g) 2 for the existence of superradiance, one obtains 
an allowed range of permissible values of N. Requesting an upper bound much larger than 
the lower bound leads to the requirement T <C M(nd) 2 /kg ~ (/t/MHz) 2 (<i/yum) 2 x 8mK in 
the case of 87 Rb. This leaves a comfortable range of the Rb temperatures for resonators with 
a quality of order 1000, corresponding to k ~ 1 MHz. Using temperatures close to the lower 
bound considered above will have the advantage that one may even record many runs with 
the same atom positions, and thus separate quantum fluctuations from the superradiant 
dynamics (e.g. waiting time distributions after initial excitation) from fluctuations due to 
changing coupling constants. As an example, for k ~ 10 6 Hz, g ~ 400 Hz, we have the 
requirement N 10 7 . At a temperature of 1 /j,K, the thermal speed of the Rb atoms is of 



the order 1 cm/s, and for d = 10 /im, we need N 3> 20. 
We now proceed with the analysis by rewriting (JS} as 

£|>] = £ [p] + eA[p] (12) 

£o[p] = ^j{[J-,pJ+] + h.c) (13) 

£ iH = ^7 E + + ^ff?) (^-V?) + ^- c -) ( 14 ) 

Yj E + *3) (V-> ^ + , (15) 



where in the last step we have neglected the terms of order 5g 2 . We have introduced a book 
keeping parameter e, assuming that C\ is small compared to Cq. In the end we will set 
e = l. We also expand p in terms of e, 

p = p + tpi + e 2 p2 • • • • (16) 
To order e°, we retrieve the original Lindblad master equation 

^ = AM (17) 
with initial condition po(0) = p(0). To first order, e , we get 

^ = £oM + £iM, Pi(0) = 0. (18) 
The formal solution of (|18|) is given by 

Pl (r)= [ T dr'e^-^de^'lpiO)}. (19) 
</o 

It describes the usual situation of propagating the initial density matrix with the "free" 
propagator exp(£ ''" / ) corresponding to the homogeneous case up to a time t', then have 
the perturbation C\ act at time r', and then continue with the "free" propagation till the 
end of the time interval. The free propagator has been well studied, and in particular very 
precise semi-classical expressions exist, both for the propagation of the probabilities and 
for the coherences — but only for j = j max relevant for the homogeneous case 23j, I25N29]. 
Here, we have the added complication that C\ does not conserve j, as C\ is by definition 
not symmetric under exchange of atoms. We therefore need a more general expression for 



the free propagator for the time interval [r', r]. We will derive here a free propagator valid 
in the entire 2 N dimensional Hilbert space, represented in all SU(2) irrep components, by 
generalizing the method in 29| that allows one to connect the propagator for coherences 
to the one for probabilities. Secondly, we will obtain explicit expressions for the matrix 
elements of L\ in and between different SU(2) irreps. Taken together, this will allow the 
construction of the full propagator according to Eq. (fl8|) . 

The irreducible representations of SU(2) can be constructed by adding one spin-1/2 after 
another. When adding a new spin, the available values of total angular momentum j, can 
either increase or decrease by 1/2. Therefore, for N > 2, there are many ways to get to 
a particular value j for given number of atoms N. These different ways lead to different 
degenerate irreps for the same j. A quantum state must accordingly be labeled not only 
by the total angular momentum j and m (m — —j, . . . , j), but also by additional quantum 
numbers a, which distinguish different irreps with the same j. We may consider a as the 
label of a path on the lattice of allowed (j, N) combinations that leads to the value of j at 
hand. We thus have states \a; jm) that are eigenstates of J 2 and J z with eigenvalues j(j + 1) 
and m, respectively, which we can use to represent C%. It will turn out that to first order 
perturbation theory, one needs only the irreps with j = j max-, 3 max — 1 and j max — 2. This is 
a consequence of the fact that C\ only contains tensor operators of rank 1 and 2. 



III. SUPERRADIANCE FOR HOMOGENEOUS COUPLINGS IN THE ENTIRE 
HILBERT SPACE 

It can be shown that the superradiant master equation ([8]) with homogeneous couplings 
conserves both j and a, i.e. 

(/3;fm'\J±\a;jm} = S fj 8 a ^S mlm±1 d ± (j,m) (20) 

with d±(j,m) = a/j(j + 1) — m(m ± 1). Conservation of j has been known from the early 
days of superradiance and allowed a formulation of the problem in the j = j max irrep. 
Conservation of a can be shown by complete induction. Below we will consider the situation 
of two sub-ensembles of atoms, each of which is homogeneous on its own. In that case we will 
see that conservation of a also follows from selection rules encoded in the angular momentum 
algebra. 
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It is worthwhile to express the matrix elements of p as 

Pm(j,f,k,a,a',r) = (a;jm + k\p(r)\a';j r m — k) , (21) 

where m is a "center of mass" quantum number, m = (mi + 1712) /2, and k = (mi — m2)/2 
for a matrix element (a;j, mi\p\a';j', 777.2). A value different from zero, j 7^ j', or a ^ a' 
signifies a coherence. The quantum numbers m and k are simultaneously integer or half- 
integer. In order to render the notation less cumbersome, we denote the set of quantum 
numbers (j,f, k) collectively as x. We can then write Eq. f fTT|) as 

^-p m (x, a, a', t) = a m+ i(x)p m+ i(x, a, a', r) - b m (x)p m (x, a, a', r) , (22) 

and see that probabilities and coherences do not mix under time evolution, nor do coherences 
defined through different combinations of j,j',a,a', and k. The coefficients a m and b m are 
independent of a, a', and are defined as 

a m (x) = a m (j,j',k) = jd-(j,m + k)d-(j',m - k) 

b m (x) = b m (j,j',k) = -^j (d + (j,m + k - l)d_(j,m + k) + d-(j',m- k)d + (j',m- k - 1)) . 

The master equation can be solved with the help of the propagator D = exp(C r), which is 
non-diagonal now only in the quantum numbers m. We can thus write 

n t (x) 

p m (x,a,a',r) = ^ D mn(x, r)p n (x, a, a', 0) . (23) 

where we have used that the propagator D does not depend on a, a'. The sum over n 
runs over all values from n^x) = max(— j — k, —f + k) to n^(x) = min(j — k,j' + k) in 
steps of 1, and m is restricted to the same interval. These bounds follow from requesting 
—j<m + k<j and —j' < m — k < j' . The initial condition is D mn (x,0) = 5 mn . The 
equation can be solved exactly by Laplace transformation, 



oc 



D mn (x,z)= e ZT D mn (x,T)dT , (24) 
Jo 

which leads to the recursion relation 

(z + b m (x))D mn(x, z) — a m _(_i(x)Z) m _|_i n (x, z) + S mn . (25) 

The solution 

l=m 
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valid for all m < n, has poles at z — — 6/(x). It remains to perform the inverse Laplace 
transform, 

■y rZQ+ioo 

D mn (x,r) = — dze ZT D mn (x,z) , (27) 

^1 J z -ioo 

where zq G K. needs to be chosen larger than the real part of all poles. Unfortunately, for 
large N, there are many poles, which makes the inverse Laplace transform cumbersome. 



29) 



However, we can relate the propagator for coherences to the propagator of probabilities 
by noting that 

k(x) = bl (j,j,0) - + 1) - f(f + 1) +2k 2 ) /(2J) (28) 

ai(x) = ai(j,j,0)qi(x) , (29) 
f , y/(j-k-l + l)(j + k + l)(f + k - I + ijg - k + f) 

qi{x) = + • (30) 

Once we shift the integration variable z in Eq. (|27|) by — (J (J + 1) — + 1) + 2k 2 ) /(2J) 
and define 

n 

Q(x,m,n) = J J q t (x) , (31) 

l=m+l 

we find 

D mn (x,r) = e^ +1 )-^' +1 ^^ , (32) 

where D mn (j, j, 0, r) is the propagator of p robabilities. A precise semiclassical approximation 
for D mn (j,j, 0, r) can be found in [25|, |29|. We thus have through Eq. f[3"2"j) immediate access 
to a propagator in the entire exponentially large Hilbert space. 

The independence of D mn (j,j',k,r) on a, a' implies that coherences p m (x, a,a',r) with 
k = between irreps with different a, a', but the same j, decay just as slowly as the 
probabilities within any irrep with the same j. This constitutes another example of slowly 



decohering Schrodinger cat states, with which the theory of superradiance is rich [30|] . The 
physical reason for this is that the Lindblad operator £o[p] is totally symmetric under 
permutation of atoms, and therefore cannot make transitions between irreps with different 
values of a. Previously known examples of slowly decohering Schrodinger cat states in 
superradiance include superpositions of angular momentum coherent states in the j = j max 
irrep that are symmetric with respect to the equator J z = 0. They also decohere just as 
slowly as the corresponding probabilities. Furthermore, there is a well known decoherence 
free subspace (DFS) of dimension (^ 2 ) ~ 2 N /y/N that contains the states \a;j, —j) of all 
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irreps (oe,j) in which all superpositions, as macroscopic or entangled as they may be, are 
decoherence free under this collective decoherence process 31]. 



and 



IV. TWO SUBSYSTEMS 

The superradiance process depends on all individual coupling constants. To simplify the 
problem at hand, let us look at only two sub-ensembles with Aq and N 2 atoms (Ni+N 2 = N), 
and coupling constants g\ and g 2 , respectively. While this exact distribution of coupling 
constants may not be very realistic, one may hope to develop an idea of the effects of 
inhomogeneity in this situation that may still be qualitatively correct for more complicated 
distributions of couplings, and to find a solution that is still analytically tractable. For an 
arbitrary distribution of coupling constants one may choose to adjust N\, g\ and g 2 to the 
first three nontrivial moments of the distribution. The results are 

(Sg 3 ) - (5g 2 ) (g) ± ^4 (5g 2 ) 3 + (Sg 3 ) 2 - 6 (Sg 2 ) (Sg 3 ) (g) + 9 (Sg 2 ) 2 (g) 2 
9l ' 2 ~ 2W) ' (33) 

Ni -n 2 = n , m^mm m 

Here, (Sg 2 ) = Yli=i(9i ~~ (d)) 2 /N, etc. We assume that iV is big enough such that rounding 
to the next integer or half integer value for N\ does not lead to a significant change. The 
advantage of working with just two different couplings is that explicit analytical formulas 
for Ci can be obtained. We define the eigenstates \ji,m) of jW.jW and Jg for i — 1,2. 
The basis functions \{jiih)3 m ) which are eigenstates of J^.J^, J^.J^, J.J, and J z , will 
be used to represent the state of the whole system. The symmetry label a corresponds then 
to the pseudo- angular momentum quantum numbers ji, j 2 for the two sub-ensembles. 

For given couplings g\ = g + 5g and g 2 = g — Sg, we have 5g = (gi — gi)/2. We also 
introduce the ladder operators J+\ J_\ J+ \ and in the two subsystems, defined as 

Ni N 
i=l i=JVi+l 

C\ in Eq. (f!5l) then becomes 

CM = f {2(J W pJ^ - Ji 2) pjf) - P(jj 1} J W - 4 2) J {2) ) - (J?J W - Jf Ji 2) )p)(36) 
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We then need to calculate the matrix elements of </_,</_, and J_ , J± J_ in basis 
states |(ji, 32)3™)- The results can be found in the Appendix. All these matrix elements 
are real. The expressions for Ji 2) and j| 2) / 2) are obtained from those for J_ and J_ 
simply by exchanging j\ -h- j 2 , which evidently leaves the expressions unchanged if j\ = j 2 . 
Equation ( Tl5l) then tells us that £i[p\ = for j\ = j 2 , for any p that is symmetric under the 
exchange of the two subsystems, 1 — > 2. This means that to first order perturbation theory 
the effect of the inhomogeneity in the couplings vanishes exactly if the two sub-ensembles 
contain the same number of atoms and the initial state is in the j = j max irrep. 

A. Full propagator 

Let us consider the case of an initially pure state, p(0) = \i/j(0))(ip(0)\ with support only 
in the irrep with j = jmax- Since we start off with a totally symmetric state, the initial state 
is also totally symmetric under permutations in each sub-ensemble, and we have therefore 
initially j% = Ax/2, j 2 = A 2 /2. To zeroth order we find 

Po{r) jm ,j'm' = dj,j'$j,j max Yl D rn±mL n U,J, 7, > T W n+ ^mL (°)^„_in^ (°) > ( 37 ) 



with V'm(O) = ((ji, j2)jmax m \4'(fy) ■ The first order correction to the density matrix can be 
written as 

Pl(r)jm,j' m ' = 5g) Pjmj>m';srr>l/>*r-r> (Q)^„ r -r> (0) , (38) 
s,r,r' 

with a propagator P defined as 

f T x ^ m — m' , . r-r' , , 

-Ljmj'm';srr' / / m+m' r l.7i .7 j ~ ? 7" 7" j r+r' „ ( Jmax ; Jmax i ~ ) 7" JOT 

JO 2 * 2 * 



L 2 . (39) 

Jmax i Jmax f 



The sums over s, r, r' are restricted such that —jmax < s + (r — r')/2, s — (r — r')/2 < j ma x- 
We have introduced the representation of Ci in the SU (2) irrep states, 

L = T (2(jm\J W \lr)(j'm'\J W \l'r')* 

I T V V 1 Jmax jj 

-^M^Vf J^b'W) - <^<W0'm| J W |/r) - ( ?1 o j 2 )) (40) 
With these expressions we can now evaluate any time dependent expectation value. 
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V. RESULTS 



A. Short time behavior of population inversion 

The population inversion is given by J z = | Y^iLi a z To zer °th order we have 

Jmax 

(Jzfi{r)) = tr (J t po(T)) = ^ m Po{ T )j max m,j max m 

TTl = Jmax 

= ^ mD mn (j max , j max , , t ) | ^» (0 ) | 2 , (41) 

mn 

and to first order 

(r)^ +r ^(0)V s _^i(0). (42) 



j,m,s,r,r' 



j m j' m! 

The propagator involves the matrix elements L now only for j' = j, m! = m, 

I r V r' 

and I = V = jmax- Since j < j max , the summation over j in ( 142]) is restricted to j = 
jmax, jmax — ^,jmax — 2 (see the remarks after Eq. ( 1A19I) ). Moreover, it is possible to obtain 
closed, compact expressions for these matrix elements, at least in the case of r = r', which is 
relevant if the initial state is a Dicke state ji^jmaxfn)- To this end it is useful to consider 
the three cases j = j m ax,jmax ~ 1 and j = j max - 2 separately. 

1. For j = j max , matrix elements of and of contribute, and we find 



Jmax H Jmax fl 2jj — j 



$n,r—l(jmax 4" j m ax ^1 ^1 



Jmax \Jmax + 1/2) 

SnAjmax + j max + n - n 2 ) . (43) 



Imax T Jmax f 

,■2 i „ 2 



2. For j = j max — 1, only J_ contributes, due to the Kronecker-deltas that come with 



jl(jmax Jl)( 1 + jmax ~\~ ^}\jmax + ^) / ■ ■ j 

{jmax + l/2)j max {2j 

max x I 



J { + ] J m . We obtain 

jmax 1.TI jmax \ TL 
Lj = A0 n r—1 

Jmax T Jmax f 

= , (44) 

where in the last step we have used ji = j max —ji such that the prefactor ji(j ma x—ji) = 
J1J2 is symmetric under j 1 -h- j 2 . 
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3. For j = jmax — 2, we find immediately that this contribution vanishes, as cannot 
change j m ax by more than one unit, and for the term the prefactor is again 

zero due to the Kronecker-deltas. 

The only contribution to (J Zj i) stems therefore from the original fully symmetric irrep 
with j = j max . Experimentally, the most relevant situation is a fully excited initial state 
Kiij 32) jmax, jmax) ', which can be achieved e.g. by optical pumping with external laser light 
that acts on all atoms in the same way. In this case, one obtains the explicit result for the 
first order correction, 

(Jz,i(t)) = 25g——^-. — ^ x j ^mD mn (j mm ,j max ,0,r-r') 

Jmax{Jmax + V^J JO nm 



(jmax + jmax n n ) ^n+ljmoi (jmax , jmax , , T ) 
-(jmax + jmax + n ~ ^) D n,j max (jmax, jmax, 0, t') (It . (45) 

Equation ( 1451) is one of the main results of this paper. It shows once more that (J Z) i(r)) 
vanishes for ji = j 2 = jmax I '2 . Moreover, all dependence on N\ is in the prefactor 2j\—j max = 
Ni — N/2, such that for given N it is sufficient to calculate (J z ,i(t)) for any Ni ^ N/2, and 
then rescale accordingly. Figure [3] shows (J x i(t)) calculated from Eq. (H2|) and from exact 
diagonalization of the full propagator for a small number of atoms and for very small 8g, 
such that (J z ,i(t)) / '5g ~ ((Jz(t))\8§ — (Jz(T))Ug=o)/Sg. The agreement is perfect. We see 
that (J Z) i(r)) vanishes at r = 0, and then decreases as function of time, before increasing 
again. This implies a more rapid initial decay of the population inversion for 5g > than for 
the homogeneous case. As the loss of atomic excitation goes along with an increase of the 
photon number in the cavity, superradiance is accelerated by the inhomogeneity. Note that 
5g > together with j\ > j 2 means that there are more atoms with the larger of the two 
coupling constants, such that one indeed expects an acceleration of superradiance compared 
to the homogeneous case. For ji < j 2 , (J z ,i(t)) changes sign, and superradiance slows down. 

Figure [3] also shows that for sufficiently large times, (J z ,i) vanishes. This can be under- 
stood from Eq. (1431) . as for large r the two free propagators under the integral cannot be 
simultaneously substantially different from zero for the combinations of indices that appear. 
As a consequence, the effect of inhomogeneity on the average value of (J z (t)) for r ^> 1 is 
at most a quadratic function of Sg, even for j\ ^ j 2 (see also next subsection). This means 



16 




FIG. 3: (Color online) First order correction (J Zj i(r))/Sg to the population inversion (J z (t)) as 
function of dimensionless time r. Green dots are obtained from the perturbation calculation, 
Eq. ()42p : blue continuous line: exact result from numerical diagonalization of full propagator for 
a small number of atoms. Curves from top to bottom correspond to (N\,N2) = (2, 1), (3,2), (3,1), 
and (4,1), respectively. 



that there is a finite time after which first order perturbation theory becomes inadequate, 
and higher order corrections will dominate. 

The time-independent factors in Eq. (14"5]) scale proportional to for j max ^> 1. How- 
ever, the integration over r brings down a factor l/(j ma x + 1/2) (see Eq. (1321) . and Eq. (4.56) 



m 



3l|). The free propagators are of order 1, such that (J z ,i(r)) / j max scales as Sgj^^ as it 
should, and first order perturbation theory remains meaningful for large j ma x- This means 
that it is enough to calculate (J z ,i(t)) / j ma x for moderate values of j ma x, as it will saturate 
as function of j ma x- 



B. Incomplete Relaxation 

Another important consequence of inhomogeneous coupling constants that can be ob- 
served even for a very small number of atoms is incomplete relaxation (neglecting spon- 



ul 



y SU(2)- 



321 ] contain- 



taneous emission — see the remarks in sec JII Aj) . It is well known that for 
symmetric superradiance there is a large decoherence free subspace (DFS) j^l, 
ing ( N ,^) dark states for N even (of ((j V _ 1 w 2 ) f° r N odd) . These states, defined through 
J-\ip) = 0, can trap the dynamics, in the sense that if such a state is reached, the superradi- 
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ant dynamics is switched off and further evolution is only possible through competing mech- 
anisms neglected so far. The simplest example is given for N = 2 with two DFS states. If 
we denote the two hyperfine states involved as |0) = \F, mp = 1, 0) and |1) = \F, rrip — 2, 0), 
the DFS states are the ground state |00) and the "singlet" state (|01) - \10))/y/2. In the 
singlet state both atoms together contain one photon, but destructive interference prevents 
the transfer of the photon from the atoms to the cavity (from where it would escape). A 
way of reaching that state is to start in an initial state |01), which in half the cases will 



emit a photon, but in the other half get trapped in the singlet state |33|. This constitutes 
a simple way of preparing an entangled state through a decoherence (and even dissipation) 
mechanism: If no photon leaves within a time given by max( 1/g, 1/ k), the system is with 
high probability in the singlet state. 

More generally, for perfect SU(2) symmetry, the DFS states are the (typically highly 
degenerate) ground states of all the SU (2) irreps with j — 0,1, ... , N/2 (assuming N even). 
If the SU{2) symmetry is broken, the DFS does not disappear, but rather gets rotated 
in Hilbert space. For example, for N = 2, and real coupling constants g\ and g%, the 



singlet is replaced by a state 5^2 |0)i| 1)2 — t7i 1 1) 1 1 0) 2/ ygT+ <?2> which is still annihilated 
by gi<j^ + gi<j_, the new collective Lindblad-operator. This parametric dependence of 
;he DFS on a system parameter is at the basis of "decoherence-enhanced measurements" 
34|, |35|, which allow precision measurements with Heisenberg-limited sensitivity while using 
initial product states. 

With perfect SU(2) symmetry, an initially fully excited state with j = N/2 = j max 
remains in that irrep and relaxes to the ground state, without ever reaching a nontrivial 
DFS state, which only exists for j 7^ jmax- However, when the SU{2) symmetry is broken, j 
is no longer conserved, and nontrivial DFS states can be reached, resulting in the trapping of 
the population. Since the first order correction to (J z (j)) vanishes for large r, the trapping 
effect is beyond reach of the perturbation theory developed above. We therefore resort to a 
numerical approach by simulating the stochastic Schrodinger equation (SSE) that unravels 
Eq. ( fTTT) . The SSE is given by 

dij){t) = D 1 (ip(t))dt + D 2 (ip(t))dW(t), (46) 
D x (V0 = 7 {2{J-hJ- ~ J+J- ~ (J-)J) i> (47) 
D 2 ty) = v / 27(^-(^)^)^, (48) 



18 



dW(t) is a Wiener process with average zero and variance dt, and (J_)^ = 36]. 
When averaging over a large number of realizations of the stochastic process one obtains a 
numerically exact solution of the master equation. In principle this can be done even for 
arbitrary coupling constants, but we stay with the situation of two different subsystems. A 
drawback of the numerical approach is that it is limited to a small number of atoms. 

Figure H] illustrates the incomplete relaxation in the inhomogeneous case for typical real- 
izations of the stochastic process from Eq. (j4"6]) . For Sg = 0, superradiance proceeds down 
to the ground state, whereas for Sg = 0.5 the trajectories can get trapped at a random 
finite value (J z (t — > oo)). The statistical analysis of these final values (we take r = 5, as 
the trajectories typically have converged at that time) leads to the histograms shown in the 
same figure, obtained from 10 5 runs of the SSE, Aq = N 2 = 2, and Ni = N 2 = 3. For Sg = 0, 
the histogram is a 5-peak at {J z {T~))/j = —1, as the superradiant relaxation proceeds to the 
total ground state. For Sg — 1, a 5-peak at (J z ( T ))/j = arises. This is due to the fact that 
in this case the second set of atoms has coupling constants zero. Therefore, these atoms 
remain excited, whereas the atoms in the first set decay to the ground state, such that in 
the end half of the excitation remains in the system, resulting in (J z ) = 0. In general, for 
Aq ^ N 2 the final value for Sg = 1 is given by (J z (t ->> oo))/j = (N 2 - iVi)/(JV a + Aq). For 
intermediate values of Sg, the histogram rapidly broadens and shifts to larger values of (J z ) 
with increasing Sg. 

Figure [5] shows the average value of (Jz{t))/j for large r (r = 5) as a function of Sg 
obtained from these 10 5 realizations and the standard deviations (as errorbars) for A = 4, 
Aq = N 2 = 2, and for N = 6, Ax = N 2 = 3. Initially we see quadratic behavior, which 
arises from the vanishing of (J Zj i(t)) for Aq = N 2 or r 3> 1. 



VI. CONCLUSIONS 



We have presented a thorough analysis of the effect of superradiance of cold atoms coupled 
to a superconducting on-chip LC-resonator. Under realistic conditions we demonstrated a 
parameter regime in which superradiance should be observable. We have analysed the effect 
of inhomogeneous couplings on the superradiance process by perturbation theory in the 
inhomogeneity, and numerical simulations. By dividing the sample into two subensembles, 
with different atom numbers and coupling constants, we can model the inhomogeneous 
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FIG. 4: (Color online) (a) Ten different realizations of the stochastic Schrodinger equation (|46p for 
5g = (black) and Sg = 0.5 (red/grey) for N = 12, N x = N 2 = 6. Plotted is (i/j(t)\ J z \ip{r))/j 
for randomly evolving states \tp(r)). (b) Histogram of (J z ( T ))/j for large r (r = 5) for N = 6, 
Aq = A^2 = 3 from 10 5 realizations of the SSE Eq. (|46|) for different values of 5g (5g = 0.1 to 
5g = 0.9 in steps of 0.1 from left to right). The J-peaks at {J z { T ))/j = — 1 and corresponding to 
5g = and bg = 1 are not shown. 




FIG. 5: (Color online) Full lines: Average values of (J z ) / j (over 10 5 realizations of the SSE Eq. (|46p ) 
and standard deviations (dashed lines) of (J z (_ T ))/j for large r (r = 5) as a function of 5^. Black 
circles: N\ = N 2 = 2; red/grey squares: N\ = N 2 = 3. Full and dashed lines are guides to the eye. 

coupling. Our results show that superradiance may be accelerated or slowed down compared 
to the homogeneous case, depending on the distribution of coupling constants. The first 
order correction in the inhomogeneity vanishes for all observables in the case of two sub 
ensembles containing the same number of atoms. For large times, inhomogeneous coupling 
constants can lead to population trapping in random decoherence free states, and we have 
provided numerical results for the population inversion starting from a initially fully excited 
state. The first order correction to the final average population inversion vanishes in this 
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case, such that the change of that quantity is at most quadratic for small inhomogeneity. 
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their computers. 

Appendix A: Evaluation of matrix elements of angular momentum operators 

We use two different ways of calculating matrix elements of angular momentum operators 
in the joint basis j2)j'm): i.) Decoupling the basis states into single angular momentum 
basis states using Clebsch-Gordan coefficients and ii.) using the Wigner-Eckart theorem. 

1. Decoupling into single angular momentum basis states 

A straight forward way of obtaining the matrix elements of J± and J + J_ is to decouple 

jl ... 
\ Jl J2 J 

\{ji,j2)jm)= 22 C |ji"ii)|j' 2 7n-mi), (Al) 

mi=-h m-mi m 

where the coefficients are Clebsch-Gordan coefficients. Then apply the desired operator, and 
couple the states back together using the inverse transformation, 

h+h ■ ■ i 

\jimi)\j 2 m - mi) = ^ C \{ji,h)lm) (A2) 

HJW2I mi m2 m 

Since the J± are symmetric under permutation of the atoms in subsystem i, they conserve 
ji, such that the matrix elements are all diagonal in the index a = (ji, j 2 ). We find 

((j 1 j2)jm\J i y\(j' 1 j' 2 )j'm') = J2 d-UnmjC Jl 32 3 ' 

mi=-ji m i rn + \ — m\ m + 1 

C 3 3 b 3 \,3>J>nj'J>m,m'-\ (A3) 

mi — 1 m + 1 — mi m 

jl ■ ■ 7 

((j l ,j2)jm\4 ) J^\(j' 1 ,j' 2 )j'm f ) = J2 d + {^ mi -l)d^ mi )C 31 32 ' 



mi=—ji 



mi m — mi m 



C 31 3 SjuAj'Jm,™' (A4) 



mi m — mi m 
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Another derivation which in the end gives closed analytical expressions is based on the 
Wigner-Eckart theorem. 



2. Wigner-Eckart theorem 

Let us consider momentarily a single angular momentum j (i.e. with Hilbert space di- 
mension 2j + 1 spanned by | jm) basis states which form a simultaneous eigenbasis of J 2 and 
J z ). The Wigner-Eckart theorem states that the matrix elements of an irreducible tensor 
operator T^q which transforms according to the irrep of SU(2) with j = K, i.e. like a state 
\KQ), is given by 

(jm\T KQ \fm>) = (-l) 2K (j\\T K \\f)C J ' K 3 , (A5) 

m! Q m 

where (j\ \Tk\ is a reduced matrix element that does not depend on the magnetic quantum 
numbers m, m! or Q [37]. In practice one calculates these by using the Wigner-Eckart 
theorem backwards for a simple operator Tkq whose matrix elements are known. There is 
just one scalar tensor operator that can be formed from the components of J, T 00 (J) oc J 2 = 
j(j + 1)1. Tensor operators of rank 1 (i.e. a vector) are formed by the components of J. We 
have 

T 1±1 (J) = t-L j ±> T 1Q {3) = J e . (A6) 

Higher order tensor operators of rank up to K = k + k' can be formed from the product of 
lower rank tensors Rk q , S^y 

„ k k' K 

TKQ(Rk, Sk') = / jRkqSk'q'C . (A7) 

q, q ' q q' Q 

One particular example is a tensor formed by the Cartesian product of the components of 
the vector operator Ti 9 (J) (q = —1,0, 1) introduced above, which we denote as Tkq(J, J), 
and which reads 

T KQ (3,3) = ^T l9 (J)7V(J)C 1 1 K . (A8) 
q, q > q q' Q 

The Clebsch-Gordan coefficients limit the possible values of K to K = 0, 1, 2. We can invert 

this relation and obtain the reduction of a product of irreducible tensor operators into a 
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sum of irreducible tensor operators, 



2 11 2 



T lq (3)T lq ,(3) = J2 C T lq+q ,(J,J) . (A9) 

1=0 q q' q + q' 

With the help of this equation we can reduce the product J + J_ into its irreducible compo- 
nents, 

J + J_ = -2T 1+1 (J)T 1 _ 1 (J) = -2 ^-Lt 00 (J, J) + ^=T 10 (J, J) + -^=T 20 (J, J)) . (A10) 

Note that Ti i3 (J, J) = —T\ >q {3)/y/2. This relation can be shown component by component 
using Eqs. (lA8)A6j) . and the commutation relations of the angular momentum operators. 

Now consider a composite system of two (physical or pseudo-) angular momenta ji and 
j 2 . We distinguish the operators acting on subsystem i as before by a superscript, T^h- The 
reduced matrix elements of an operator which acts only subsystem 1, Txq = T^q eg) l( 2 ), 
can be related to the ones in subsystem 1 alone according to 

((ji,j2)j\\T K \\(j[,j' 2 )f) = (-l)' + 'i-*-*vW + l)(2 Jl + Dli'!./:. ./•;../../' K, j 2 )( 3i \\t£%[}6 J2 ,, , 

(All) 

where the symbol W is related to Wigner's 6j-symbol by 

\V(u.l,r.<l:r.f) = (^iy+ b + c + d { ab 6 } . (AI2) 

d c f 



see Eq. (5.9) in [37|. From ( jAlOlATTjl we obtain 

(Uuj2)jm\J { } ) \(j 1 ,j 2 )j'm') = V2(-l) j+h - h+ W(W + l)(2ji + 1) 



l|-(./,../.../../':l../2K' / 1 ' (jillT^Wibi} (A13) 
m' — 1 m 



((ji,j2)jm\J^J W \( 3l ,j 2 )fm') = 2(-l)^^ + V(2j' + l)(2ji + l) (A14) 

{i=^(j 1 ,j 1 ,j,j";0,j 2 )C y J O'iII^CJ.J)!^) 

V J m' m 



-—wUi,jiJ,f;i,h)C (J>J)lb'i) 

v 2 m' m 

+ * W (j 1 ,j 1 ,j,j , ;2,j 2 )C J 2 3 (j 1 \\T^(3,3)\\j 1 )\. 
Vo m' m J 
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It remains to calculate the reduced matrix elements. Brink and Satchler [37| (p.51ff) give for 
a single angular momentum Too(J, J) = J 2 /V3. Using this in the Wigner Eckart theorem, 
sandwiched between states (j0\ and \j'0) gives for a single angular momentum 

(j\\T (3,3)\\j) = -^=j(j + l). (A15) 

Similarly, from Tio(J) = — J z we find 



(m(j)\\j} = vM+i), ( Al 6) 

and 

(j||T 1 (J,J)||j) = - v / J ( J + l)/2. (A17) 

Finally, one can show that T 2 o(J, J) = (3 J 2 — J 2 )/\/E by applying ( 1A8I) . and the Wigner- 
Eckart theorem tells us that hence 



(j||T 2 (J, J)||j) = v/j(j + 1)(2j-1)(2j + 3)/6, (A18) 

and thus 

(Ui,j2)jm\j w \u 1 , 32 )fm') = V2(-iy + ^ + W(W + l)(2ji + + 1) 

\V(.n.J..J.j':\.jM' J 1 3 (A19) 
m' —1 m 

(Ui,h)jm\J^ ) J^\(ji,j2)fm') = 2(-l)^ 1 -^ + V(2j / + l)(2j 1 + l) (A20) 

{- \wUi,ji,j,f;0,j 2 )c 3 ° 3 ii(ii + 1) 

" m' m 

+V(ii ) ii,ij";i I i2)c J 1 J ViiO'i + i) 

m' m 

+lwUi,ji,j,f',2,j 2 )C J 2 3 v / Ji(ji + l)(2ji-l)(2j 1 + 3)}. 

m' m 

We have checked numerically up to j = 10 and all combinations of j',ji,m,m' (with j 2 = 
f — ji) that Eqs. (lA19|A20p give exactly the same results as Eqs. (IA3IIA4j) . respectively. 
Eqs. ( lA19fA20j) have the advantage of avoiding an additional sum, and of leading to explicit 
expressions for the case relevant for the calculation of (J z (r)), as we demonstrate in Results 



24 



section. Furthermore, the selection rule j G {j'—2, . . . , j'+2} is evident from these equations, 
as otherwise the Clebsch-Gordan coefficients vanish. 
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